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Abstract 

Marching surfaces is a method for isosurface extraction and approximation based on a multi¬ 
sided patch interpolation scheme. Given a 3D grid of scalar values, an underlying curve network is 
formed using second order information and cubic Hermite splines. Circular arc fitting defines the 
tangent vectors for the Hermite curves at specified isovalues. Once the boundary curve network is 
formed, a loop of curves is determined for each grid cell and then interpolated with multi-sided surface 
patches, which are G^ continuous at the joins. The data economy of the method and its continuity 
preserving properties provide an effective compression scheme, ideal for indirect volume rendering on 
mobile devices, or collaborating on the Internet, while enhancing visual fidelity. The use of multi¬ 
sided patches enables a more natural way to approximate the isosurfaces than using a fixed number of 
sides or polygons as is proposed in the literature. This assertion is supported with comparisons to the 
traditional Marching Cubes algorithm and other G^ methods. 

Categories and Subject Descriptors (according to ACM CCS): G.1.2 [Numerical Analysis]: Numerical 
Analysis—Approximation of surfaces and contours L3.5 [Computer Graphics]: Computer Graphics— 
Curve, surface, solid, and object representations L4.8 [Image Processing And Computer Vision]: Scene 
Analysis—Surface fitting 


1. Introduction 

Isosurface extraction, approximation and display is 
fundamental in visualization for volume data, e.g., 
in scientific visualization. Marching cubes [LC87] has 
been the de facto standard for indirect volume ren¬ 
dering since its inception in 1987; major scientific soft¬ 
ware still uses it despite disadvantages such as serrated 
edges, angularities, linear precision, and particularly 
large polygon count. 

Our goal is to create an isosurface representation 
and rendering method that is visually appealing, yet 
accurate, useful for mobile computing and collabora¬ 
tion, and ultimately able to leverage recent parallel 
hardware capabilities in a parallel version that han¬ 
dles large datasets in accordance with [Joh04]. This 
paper introduces the method and demonstrate its cor¬ 
rectness and compressibility. We selected high-order 
surfaces that provide smoother reconstruction and ex¬ 
hibit fewer rendering artifacts than existing methods. 

The first step is to reduce the 3D gridded dataset 


to a a representation on the three independent ax¬ 
ial planes. That allows us, to reformulate the problem 
of finding isosurfaces to that of finding isolines. Next, 
we deduce the connectivity of the isolines. We simplify 
this process by allowing only cell-by-cell cube-to-plane 
intersections bounded by only four different types of 
high-order surfaces, i.e. with three, four, five and six 
sides. By using neighboring cell information we pro¬ 
vide a more accurate reconstruction that avoids some 
of the well-known ambiguity cases of marching cubes 
[NH91] and allows us to have curve network with G^ 
continuity. Geometric continuity is a major advantage 
of using high-order composite surface reconstruction. 
We provide a set of conditions that enable multi-sided 
surface patches to join with G^ continuity at the cross 
boundaries. 

Related work. This paper builds on a significant 
amount of previous research in indirect volume render¬ 
ing and surface approximation methods using high- 
order interpolation surfaces. These methods can be 
classified according to the type of surfaces in which 
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they represent the isosurface or implicit function of 
interest. 

The seminal work of [GN89] on three-dimensional 
finite element simulations and other coarse volumes 
uses bi-cubic polynomials in the form of Ferguson 
patches. Their approach yields surfaces and surface 
normals that approximate high-order threshold sur¬ 
faces directly at the corner values of the cells of in¬ 
terest. One drawback is the noticeable discrepancies 
between the surface shape and light source shading 
near the bi-cubic boundaries. The authors call this 
issue ’’result smoothing”, which comes from the fact 
that Ferguson patch is not the ideal type of surface 
for this purpose. 

The work of [SZ07] describes a scheme based on cu¬ 
bic splines on type-6 tetrahedral partitions of volu¬ 
metric grids. The authors solve an optimization prob¬ 
lem to find the appropriate coefficient combinations 
that generate a boundary description of the tetrahe¬ 
dron; operators are able to provide an error bound and 
construct the tri-variate function that will serve as the 
basis for volume reconstruction. This technique falls 
into the class of quasi-interpolation methods for tri¬ 
variate splines. Related methods are [KZ08], [KKG09], 
and [MKRG12]; they leverage the relatively low to¬ 
tal degree of the splines to construct a hardware ac¬ 
celerated implementation for real time reconstruction 
of isosurfaces. Tri-linear interpolation is a common 
choice of high-order approximation methods; unfortu¬ 
nately, smoothness of adjacent splines is not a trivial 
task. Furthermore, these schemes require high-order 
derivatives that have to be approximated. 

Trimmed surfaces of triangular rational cubic Bezier 
patches is the basis of [The02], the author introduces 
it as a direct extension of marching cubes provid¬ 
ing a topologically exact approximation of contours. 
It leverages continuity at the expense of a global 
re-parametrization step at the end of its method. Al¬ 
though trimmed representations allows the use of a 
single type of patch, it adds unnecessary complexity 
compared to our method that naturally chooses the 
best type of patch cell-by-cell. 

Other related methods that provide continu¬ 
ity are [FH12], [LM08], [YNP97] and [Man92]; these 
methods achieve continuity of the adjoining control 
points by an adjustment to the boundary curve that 
ensure continuous tangent plane everywhere. While we 
make use of the same tangent plane requirement, we 
introduce the concept of tangential ribbon which puts 
no restriction to the underlying boundary curve net¬ 
work. This is the as well last step of our method, but 
with the difference that the computation is local with 
respect to every cell. 

The need of a lightweight volume rendering for re¬ 


stricted environments such as the web browser or mo¬ 
bile devices is becoming evident in recent publications 
such as [RA12], [JKD*12], and [GAR13]. The standard 
choice for indirect volume rendering is still marching 
cubes; in this paper we advocate for the high reso¬ 
lution representation here presented as the method of 
choice due to its compact representation and the econ¬ 
omy of storing and transmitting complex surfaces that 
are broken down into multi-sided composite surfaces 
with G^ continuity. 

2. Boundary Curve Network 

The first step towards marching surfaces is to con¬ 
struct an underlying boundary curve network per each 
axial plane; namely, xy, yz and xz. We extract a set of 
curves that represent isolines at a function values c, 
called isovalue. This set of points in the domain map 
to {a: e M : f{x) — c} [GJ89], as illustrated in figure 1. 
We find an approximation of x : f{x) = 0 based on cu¬ 
bic Hermite splines (1); appropriate choose of tangent 
vectors mi and m 2 (1) assures global G^ continuity 
across a connected sequence of curves and a high or¬ 
der parametric interpolation [Far02], see figure 2. 
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Figure 1: Axial plane with funetion values at eaeh 
grid point Fij = c, and an isoline representing a 
boundary curve. 

p{t) = (2t^ — 3t^ + l)po + (t^ — 2t^ + t)mo (1) 
+ (—2t + 3t )pi + (t — t )mi 

Each piece of the isoline is a third-degree polynomial 
specified by two control points and the first deriva¬ 
tives at the start and end points, over the unit interval 
(0,1); given po at t = 0 and pi at t = 1 with t G [0,1]. 
Gontrol points are placed at lines with zero crossing 
based on a linear proportion of f{x)ij values at grid 
points. Tangents are computed by fitting a quadratic 
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Figure 2: A cubic Hermite curve is defined by two 
points and two tangent vectors. Appropriate choose of 
tangents between Hermie splines assures continu¬ 
ity. 

curve trough three points at the time over neighboring 
cells, as shown in figure 3. 


p. 




Figure 3: Choosing three consecutive points to make 
use of second order information. The tangent vector rh 
at the middle point is computed by fitting a quadratic 
curve. 


Parabolas are the first curve of choice [The02] to 
derive tangents, but considering that isolines should 
stay as close as possible, we suggest that circular arcs 
provide maximum entropy (min max of curvature) as 
derived in [Str07]. It is shown that the best approxi¬ 
mation after a straight line is a circular arc; which is 
the shortest curve with minimal area (2). 


minimize 

P{u) 

subject to 



u{x) dx 


( 2 ) 

A 


Figure 4: Connecting 3 consecutive points to form a 
triangle in order to find the radius R of the enclosing 
circle, a is the length between Pi and P 2 , b between P 2 
and Pi, and c between Pi and P 3 . 


s — — (a + 6 + c) 

(3) 

K = yj {s{s - a){s - b){s - c)) 

(4) 

II 

(5) 

Using the parametric equation of the circle it is pos¬ 
sible to find the tangent vector m of the middle point. 

with angle 0 G [0, tt] and vectors a = Pi 
b = P3 — P2, see figure 5. 

— P2 , and 


(6) 

.—RsiniO). 

"^=Ucos(0)' 

(7) 


The construction of a C^ continuous set of Hermite 
curves using three consecutive points requires precise 
definitions of tangents at the first and last point. Every 
point in the isoline is composed by the tuple {tk,Pk) 
for k = 1,... ,n. For k = 1 tangent is defined by (8) 
and for /c = n by (10); tangents of the middle points, 
i.e. /c = 2,..., n — 1 are defined at (9). An special case 
has to be defined when the isoline does not form a 
loop, here we provide tangents definitions for /c = 1 
and for /c = n (12). 


A circular arc can be fitted by computing the cir- 
cumradius R (5) of the enclosing circle that is form 
by connecting the three consecutive points of inter¬ 
est with lengths a, b, and c that form a triangle; then 
the calculation of its semi-perimeter s (3) and area 
K (4) according to Heron’s Formula and [Joh29] as 
illustrated in figure 4. 


mo(Pn, Fi, P2) 

(8) 

mk(Pk-i, Pk, Pk+i) 

( 9 ) 

'kTin (Pn —1? Pni Pi) 

(10) 

mo * (Pi, P2, P3) 

(11) 

TTln * (Pn —2 5 Pn —15 Pn) 

(12) 
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Figure 5: Circular arc approximation given 3 points. 
Vectors a and b are used to find the angle 0 between 
them, to ultimately define the tangent vector rh at the 
middle point P 2 . 



Figure 6: Boundary curve network of implicit func¬ 
tion — 1 = 0 over [-3, -3, -3] x [3, 3 ,3]. 

Isovalue = 0.0. In blue are the boundary curves, in 
orange are the control points. Showing one isoline per 
axial plane. 

The length of the tangents can be adjusted for dif¬ 
ferent fits. When the adjacent points are nearly co- 
linear we suggest that 1/3 of the distance between 
endpoints is an effective factor, which comes from the 
definition of a cubic Bezier spline, that mimics the 
legs of its control polygon. As the tangent legs be¬ 
come right angles a factor of | sin(^) lengthens them 
and the curve straightens up, this gives optimal visual 
fidelity in most of the cases. 

Once points and tangents are defined at every zero 
crossings in all the cells of the plane, curves can be 
rendered in each axial plane. Up to this point we have 
an intermediate visualization method readily available 
formed out of boundary curves. Figure (6) and (7) 
illustrate the reconstruction of two implicit functions 
with its corresponding isovalues on R^. 


3. Multi-sided surface construction 

The resulting intersections across boundary curves 
naturally suggest the number of sides of the necessary 



Figure 7: Boundary curve network of implicit func¬ 
tion constructed of the blend of a sphere and a hyper¬ 
boloid |((x+l)^ + (y+l)^-l-(z+l)^-l) + |(xyz-3) = 0 
over [-3, -3, -3] x [3, 3 ,3]. Isovalue = 0.0. Show sev¬ 
eral isolines per axial plane. Isolines are globally 
continuous based on Cubic Hermite Splines. 

surfaces to fill in cell-by-cell. As mentioned before, our 
method selected a type of surfaces flexible enough to 
handle any number of surfaces sides. If we consider the 
original look up table of 15 cases [Che95] of cube-to- 
plane intersection, we can group all the possible cases 
by sets of three, four, five and six-sided surfaces only, 
as summarized in figure 8. 


Three-sides 









K 







Three and four-sides Three and five-sides 






Figure 8: Out of the 15 original lookup table cases of 
the Marching Cubes algorithm, we group the choice of 
higher order surfaces to four types of three, four, five 
and six sides. 


Surface construction is based on [GR05] and 
[VRSll]. They define an interpolation scheme based 
on the weighted, discrete least squares solution x = 
{x,y,z) that minimizes: ~ ^0^ + (2/ - Vi)^ + 

(z — Ziy]Wi(x,y, z). The least squares solution slews 
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toward the points with larger weights Wi{u). Weight 
selection for multi-sided patches is based on a para¬ 
metric map of control polygons on with the same 
number of sides k as the required surface, called foot¬ 
prints and denoted by Ui (figure 9). Boundary curves 
are then used to map each footprint’s point Ui into the 
resulting surface in R^, denoted by the attribute func¬ 
tion fi{ui). Weights are given as reciprocal distances 
from the boundary curves fi to the footprints Ui. The 
following equations (13) give the complete surface con¬ 
struction definition: 



of curves as shown in the three-sided patch at figure 
10. The complete reconstruction is shown in figure 11. 
The set of neighboring surfaces yield to a fair surface 
reconstruction with G° continuity. 



Figure 10: Partial reconstruction showing one three- 
sided surface. Arrows represent the tangent vectors at 
the adjoin points. 


Figure 9: Footprints for three, four, five and six-sided 
surfaces. 


F{u) = (x, y,z) = Y, fi{ui)Wi{u) (13) 

i 

Wi{u) = Wi{u)/Y,Wk(.u) 

k 

The previous surface definition together with a 
number of unique points per cell allows us to start 
filling in the regions that will constitute the partial iso¬ 
surface in a composite and independent manner, from 
here the selection of the word marching to our method, 
in comparison with the industry standard marching 
cubes. 

Depending of the vertex valence, there might be 
multiple tangents directions to be picked up per ad¬ 
joining point. Selecting the correct tangent direction 
is key to get a correct reconstruction. In practice, the 
implementation could use a dataset that stores the 
tangent directions and its corresponding axis, or take 
the dot product between the normalized tangent and 
the normalized vector between start and end point, 
selecting the tangent with the smallest angle as the 
appropriate. Tangent directions already suggest an ar¬ 
rangement of points in loop form, also for practice ap¬ 
plication, we suggest an additional order selecting a 
consistent wind direction for all the surface faces. The 
resulting arrangement will be the in the form of a loop 



Figure 11: Final composite surface of a unit sphere 
with eight three-sided surfaces and G^ continuity; im¬ 
age also depicts its corresponding boundary curve net¬ 
work. 


4. continuous composite surfaces 

continuity requires continuous tangent planes be¬ 
tween neighboring surfaces, and they are accomplished 
by tangential ribbons. 

Every boundary curve has a ribbon. Ribbons are 
constructed by taking the two tangents of adjoining 
curves and using them to offset the curve. The surface 
will approach the curve with the same slope as the rib¬ 
bon if the two ribbons on a curve (one for each patch) 
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are cotangent, then the surface will be. Furthermore, 
if the ribbons are the same interpolant of the tangents 
then one just has to make the tangents collinear to 
get across boundaries. An illustration of a ribbon 
of two neighboring three-sided surfaces is in figure 12. 
Figure 13 shows a complete set of ribbons for a surface 
between two contiguous three-sided surfaces. 



Figure 12: A tangential ribbon is formed as a lofted 
boundary eurve towards the direction of the tangents 
of the the two adjoining curves. Black line indicate 
that the ribbons and the boundary curve are co-linear. 
Shaded version to show the underlying triangulation. 



Figure 13: A complete set of ribbons for two neighbor¬ 
ing multi-sided patches, ribbons in the plane xy define 
tangent planes across the common boundaries and that 
yields continuity. 


Adding ribbon information to the previous multi¬ 
surface formulation (13), requires a distance function 
Si = di(u) from the set of points Ui to the footprint 
fi.Qi are the lofted cubic Hermite splines defining the 
ribbons. Finally, the weight function is squared to sat¬ 
isfy continuity requirements, as stated in [GR05]. The 
new formulation is as follows (14): 


F{u) = ^ Ri{si,ti)Wi{uf (14) 

i 

Ri{si,ti) = (1 - Si)fi(ti) + Sigi(ti) (15) 

The following figure shows two adjacent tree-sided 
surfaces with the previous G° formulation and the just 
stated G^ formulation using tangential ribbons; a con¬ 
siderable upgrade on surface approximation can be no¬ 
ticed, see figure (14). 



Figure 14: Two neighboring three-sided patches, with 
CP and G^ continuity respectively. Low resolution sur¬ 
faces disclose the underlying triangularization. 


5. Results 

A set of implicit functions are here exposed to demon¬ 
strate the output of the marching surfaces algorithm. 
Isosurfaces constructed out of three, four, five and six- 
sided surfaces make up the hnal composite surface 
of the contour of interest. Gomparisons between the 
marching cubes and the G° and G^ representations 
are shown. 

Figures 15 shows the scalar field s{x,y,z) =x^-|- 
y^-\-z^ — 1 = 0 over a [-3, -3, -3] x [3, 3 ,3] grid and iso¬ 
value = 0. Figure 15a shows the isosurface extraction 
by the traditional marching cubes algorithm composed 
of 8 triangles. Figure 15b shows the same scalar held 
approximation using marching surfaces with 8 three- 
sided surfaces and G° continuity. Figure 15c shows the 
same isosurface using ribbons with G^ continuity. 

Conclusions. We have presented a method for 
approximating and rendering isosurfaces for implicit 
functions or scientihc datasets. Our algorithm gener¬ 
ates composite surfaces out of three, four, hve and 
six possible sides. Visual appealing is directly related 
with the degree of smoothness of the approximation, 
and we described two alternatives that yield to and 
G^ continuity across neighboring edges; both alterna¬ 
tives can be computed independently, i.e. cell-by-cell. 
The ability to store large amounts of geometry based 
on multi-sided surfaces and Hermite splines yields to 
a compact and straight forward storage method that 
is amicable with restricted environments such as web 
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Figure 15: Scalar field s{x,y,z) = — 1 = 0 over [-3, -3, -3] x [3, 3 ,3], isovalue = 0. From left to 

right: (a) Marching cubes, (h) Multi-sided contours, (c) Multi-sided G^ contours. 


browsers and mobile devices. Large datasets are also 
good candidates for applications of the method here 
presented combining parallel paradigms and out-of- 
core techniques. These are all possible avenues for fu¬ 
ture work. 
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